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ABSTRACT 


The goal of this work is to develop methodologies for synchronizing the clocks of 
neighboring nodes of an underwater acoustic network. Clock synchronization requires an 
estimate of the signaling time delay between the nodes, which is particularly challenging 
in the presence of multipath propagation through the acoustic communication channel. 
This thesis focuses on modeling the underwater acoustic communication channel, 
accounting for the multipath arrivals, and creating a set of signal processing algorithms 
for estimating the required delay times that enable clock synchronization protocols for the 
underwater acoustic network. The proposed method involves correlating the responses of 
the bidirectional channels to exploit the underlying reciprocity. Performed in two stages, 
a sequence of probe signals is first transmitted to create an ensemble, which contains 
information about the time-variability of the acoustic communication channel with 
multipath. From this ensemble, we determine its dominant time-invariant characteristic 
and use it as a reference datum for the time delay measurements. The second stage 
consists in performing time-delay estimation of two probe signals exchanged between 
nodes. The two stages are tested using simulated signal measurements, and actual signal 
measurements were performed in a fresh-water lake for the first phase only. Both 
computer simulations and experimental results show the effectiveness of the proposed 
techniques. 
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EXECUTIVE SUMMARY 


Clock synchronization among network nodes becomes necessary when data from 
multiple nodes are to be co-processed. In the context of surveillance networks, 
synchronization of the nodes’ clocks is directly related to the accuracy of the localization 
of a potential target. Clock synchronization is achieved by evaluating and then correcting 
the skew and offset of the internal clock of a single node. The skew of a clock is the 
oscillator frequency difference between two synchronizing clocks, whereas the offset is 
the time difference at a specific point in time between two clocks. 

Clock synchronization is especially important for high-latency communications 
such as those in underwater acoustic networks in which data from separate nodes take 
variable amounts of time to reach the main processing node. The challenge in 
underwater acoustic clock synchronization lies in correctly estimating the time delay 
between the transmission of a signal from a node and the reception of that signal by 
another node. For the acoustic communications problem, the propagation medium 
introduces environmental factors such as high latency, multipath, scattering, refraction, 
transmission loss, noise and non-reciprocity. These factors have significant effects on the 
channel impulse response. The individual multipath arrival amplitudes and the arrival 
times can also vary greatly depending on the surface conditions and the directivity of the 
transducers. For some environmental conditions, the first or other multipath arrivals may 
be greatly attenuated or absent altogether. 

The aim of this thesis is to examine how time-variability of the impulse response 
impacts clock-synchronization and how it can be avoided or mitigated in a spatially static 
underwater network. The Time Synchronization for High Latency (TSHL) protocol 
developed by Syed and Heidemann in [1] is used as the basis for this research. The 
protocol is adapted to account for time-variability of the channel impulse response. 

In the absence of a fit-for-all underwater statistical model, channel soundings are 
a good alternative to the lack of prediction capability. They are direct measurements of 
the time-varying impulse response of a propagation medium and replace general 



assumptions about the prevailing channel between two nodes. A series of pulses is sent 
through the channel at a set interval. The receiver analyzes the received series and 
compiles its own local channel characteristics from it. 

Our approach relaxes the assumption in TSHL of a stable propagation delay 
between the synchronizing nodes. It also relaxes the reciprocity assumption for the 
channel. Instead, we repeatedly measure the response of the channel over a period of 
time to statistically define the communication channel to better estimate the propagation 
delay and its skew. This approach could also be defined as adding channel sounding to 
the TSHL protocol. Our model keeps the two-stage approach developed by TSHL, 
correcting for the skew first followed by the offset. 

In the skew-correction stage, we analyze a series of chirp responses to identify the 
least time-varying multipath present in the channel between the two nodes. Based on the 
known chirp transmission interval, linear regression is performed on the measured 
interval of arrivals of that identified multipath detected within each chirp response. The 
slope of the resulting best-of-fit straight line represents the skew of the clock’s receiver 
node compared to the originating or reference node. 

The second stage consists of an exchange of time-stamp messages between the 
two nodes, where we assume a time-varying, one-way propagation delay between the two 
nodes. The challenge in this stage consists in determining the difference in propagation 
delays. This is accomplished by performing, at one of the nodes, a time-delay estimate 
using the chirp responses of the channel sounding of both one-way propagation paths. 

The Bellhop and VirTEX underwater environment models are used to simulate a 
simple communication channel in Matlab. Variability is introduced in the model by 
modifying the state of the surface and rotating one of the nodes between the 
transmissions of the chirps. The aim of the simulations was to prove our conceptual 
approach prior to conducting an in-water experiment and to demonstrate the effect of 
time-variability on our clock-synchronization process. The simulation verified the theory 
of our approach as it precisely determined the difference in propagation delays and clock 
offset. 


xviii 



The first stage of the algorithm was then tested in the Del Monte Lake situated on 
the NPS campus. The general concept of the experiment consisted of transmitting and 
recording a series of chirps between two acoustic modems. The characteristics (duration 
and modulation) of the chirp were modified between the transmissions of the transmitted 
chirp series to analyze the impact of the chirp selection on the protocol. The digital 
recordings were recovered and analyzed using the same algorithms exercised in the 
simulation. The second stage of the algorithm was not tested since it would require 
modifying the firmware of the acoustic modems to program the message exchange 
protocol. 

Overall, the results from the experiment demonstrate some time variability in the 
tested channel. The first arrival is always present and strong. However, the following 
arrivals are observed to be gradually less stable as they increasingly interact with the 
channel boundaries. The algorithm is able to identify various multipaths based on the 
detected arrivals and detennines the most stable one based on the correlation coefficient 
of a model fit to the time-of-arrival estimates. 

Additional experiments are required to fully prove the developed algorithms and 
select the optimal parameters, such as the probe signal characteristics. Once completely 
tested, the addition of this capability to an underwater acoustic network should allow it to 
keep the clocks of the nodes synchronized. This synchronization should improve the 
reliability of the network when deployed or used for long periods of time. 

[1] A.A Syed and J. Heidemann, “Time Synchronization for High Latency 

Acoustic Networks,” in Proc. Of IEEE Conference on Computer Communications 
(Infocom), Barcelona, Spain, Apr. 2006. 
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I. INTRODUCTION 


A. REQUIREMENT FOR CLOCK SYNCHRONIZATION 

Clock synchronization in distributed sensor networks has been a research topic for 
many years. Synchronization among network nodes becomes necessary when data from 
multiple nodes are to be co-processed. In the context of target detection and 
triangulation, the synchronization of the nodes’ clocks is directly related to the accuracy 
of the localization of a potential target. 

The majority of recent developments in the literature involve wireless radio 
networks where signal propagation is almost instantaneous (speed of light). In most of 
the literature, clock synchronization is accomplished via the exchange of time-stamped 
messages between two nodes. This is achieved by evaluating and then correcting the 
skew and offset of the internal clock of a single node. The accuracy of clock 
synchronization is directly related to the accuracy of the measured propagation time 
delays among the nodes. 

B. CLOCK SYNCHRONIZATION IN THE UNDERWATER ACOUSTIC 

CHANNEL 

Clock synchronization is especially important for high-latency communications 
such as those in underwater acoustic networks where data from separate nodes take 
variable amounts of time to reach the main processing node. The challenge in 
underwater acoustic clock synchronization lies in correctly estimating the time delay 
between the transmission of a signal from a node and the reception of that signal by 
another node. For the acoustic communications problem, the propagation medium 
introduces environmental factors such as high latency, multipath, scattering, refraction, 
transmission loss, noise, and non-reciprocity. These factors have significant effects on 
the channel response. Prior to perfonning clock synchronization, we must address the 
complex time delay estimation issues associated with the response of the channel. The 
idealized model for the channel impulse response is a sequence of impulses representing 
the multipath arrivals caused by the channel. In theory, we can estimate the time delay 


1 



between transmission and reception by identifying a reference arrival (e.g. the first arrival 
or another reliable arrival). However, the individual multipath arrival amplitudes and the 
arrival times can vary greatly depending on the surface conditions and the directivity of 
the transducers. For some environmental conditions, the first or other multipath arrivals 
may be greatly attenuated or absent altogether. This makes the identification of a 
reference arrival very difficult. 

While the previously stated difficulties associated with clock-synchronization 
protocols in an underwater network are important, the extreme bandwidth and latency 
limitations of underwater communications constrain the implementation of 
synchronization protocols that involve negotiation or arbitration. Most underwater 
acoustic networks can only achieve very low data communication rates and care must be 
taken to limit the amount of information exchanged between the nodes of the network. 
Associated with this limitation is the finite power supply at each node. To maximize the 
operational utilization of a node, the number of transmissions and their associated 
message length must be limited. In the context of clock synchronization, the 
communications exchanged between the synchronizing nodes must be minimized, while 
the length of time between the need to repeat the clock synchronization algorithm should 
be maximized. 

The aim of this thesis is to examine how time variability of the impulse response 
impacts clock-synchronization and how it can be mitigated in a spatially static 
underwater network. The Time Synchronization for High Latency (TSHL) protocol 
developed by Syed and Heidemann in [1] is used as the basis for this research. It is often 
referred to as the source of other work in the literature, like this research. 

Originally, TSHL is designated for spatially fixed underwater networks. In this 
research it is adapted to account for time variability of the channel impulse response and 
tested using computer simulations. Actual acoustic signal measurements acquired under 
controlled conditions in a fresh-water lake serve to verify the physical assumptions that 
are the basis of the protocol. 


2 



This thesis is organized as follows. The general concepts associated with clock- 
synchronization and the variability of the underwater communication channel are 
introduced in Chapter II and III. Our proposed clock-synchronization algorithms 
accounting for time-variability in the channel are presented in Chapter IV. These 
algorithms are then simulated in Chapter V and some of them tested in an actual 
underwater environment in Chapter VI. Finally, the recommendation and conclusions are 
given in Chapter VII and VIII. 
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II. CLOCK SYNCHRONIZATION 


A. THE COMPUTER CLOCK 

A computer clock is an ensemble of hardware and software components used to 
provide an accurate, stable and reliable time-of-day function to the operating system and 
its clients [2]. From a hardware point of view, most clocks are based on a mechanism 
that oscillates at a precise constant frequency. The precision of the oscillator to generate 
this frequency is defined as its tolerance, and since it is usually very low, it is expressed 
in terms of parts-per-million (ppm). How well it can maintain this set frequency is 
defined as the stability of the clock. The oscillating mechanism can take the fonn of a 
simple mechanical pendulum, or the vibrations of electrons in an atom. These two 
technologies are situated at opposite ends of the spectrum in terms of precision, stability, 
complexity and cost. In the vast majority of cases, the middle-of-the-road solution is to 
use a man-made quartz crystal as the oscillating mechanism. 

By knowing the oscillating frequency of the mechanism, a unit of time can be 
linked to a set number of oscillation periods. In the case of a computer clock, a simple 
counter will send an interrupt message incrementing the clock by one time unit following 
a pre-determined number of oscillations. Sending this interrupt resets the counter to 0 
and reinitiates the process. 

B. THE QUARTZ CRYSTAL OSCILLATOR 

1. General Description 

More than 2 billion quartz crystal oscillators are manufactured annually [2]. They 
are used in consumer electronics such as wristwatches, mobile phones and computers, 
and in equipment requiring greater precision such as test and measurement instruments. 
Quartz crystal is abundant in nature and is easy to grow in large quantities, at low cost, 
and with relatively high purity. It is the piezoelectric property of the quartz crystal which 
makes it a material of choice for use as a stable resonator in oscillators. Its high Q value 
makes it an excellent material for use as an oscillator [2]. A quartz crystal can be 

manufactured in various structures and cuts, which define its own operational 
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characteristics. One such characteristic is the resonant frequency which can be from a 
few tens of kilohertz to tens of megahertz. 

2. Quartz Crystal Instabilities 

A quartz crystal often includes impurities which lead to variations in the 
resonance frequencies. This is especially true for less expensive oscillators found in most 
electronic devices, where the frequency can vary up to 40 ppm between two clocks [3]. 
The resonating frequency of a quartz crystal is subject to many types of instabilities 
which can affect its precision and accuracy [4]. These include aging, noise, and 
frequency changes depending on temperature, acceleration, ionizing radiation, power 
supply voltage, etc. By knowing how these instabilities affect a clock, it is possible to 
correct for them using hardware or software solution. As an example, the real-time clock 
used on Benthos Teledyne acoustic modems is temperature compensated [5]. 

C. THE NEED FOR A REFERENCE TIME 

Every single device in a network will have at least one local clock, the time of 
which will be different from another unless we provide proper synchronization. 
Depending on how much the times differ, the following fundamental operations 
performed by a node in a wireless sensor network may or may not be possible [3]. 

There are a number of issues strictly related to clock management, and they are 
listed next. 

1. Data Fusion 

This is one of the principal benefits to be derived from a sensor network, where 
the goal is to observe and report on events. To coherently integrate and process event 
information from multiple sensors in a meaningful way, all data must be referenced to a 
standard time. 

2. Power Management 

One way to manage and save energy in a network is by the implementation of a 
duty cycling protocol. Instead of leaving a node operational for a continuous period of 
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time, it can be programmed to sleep and wake-up at a set time interval. While in sleep 
mode, the node saves energy by keeping only essential systems operating and expending 
minimum power. To permit communication between two nodes, both of them need to be 
awake at the same time. Here lies the importance that both nodes operate on a 
synchronized clock, such that they can simultaneously switch between operating modes. 

3. Transmission Scheduling 

Some network protocols allow a node access to a network by implementing a set 
schedule for when a specific node can transmit (i.e., the time-division multiple access 
protocol) [3], From this schedule, only one node can transmit at a time, reducing the risk 
of collisions within the network. Another advantage of such a protocol is the reduction of 
the message overhead that is usually required to control the communication flow within 
the channel, therefore increasing the available bandwidth. The main requirement for this 
kind of protocol is that every node in the network must use synchronized clocks to start 
and stop their transmissions on a common distributed schedule. 

4. Master-slave Relationship 

The above operations have a mutual requirement: the internal reference clock of 
every node in the network must be synchronized. For synchronization to exist, the 
concept of a master-slave relationship must be introduced in the network. In the case of a 
sensor network, a node is designated as master and its internal clock is used as a 
reference by all the other nodes, designated as slaves. To avoid having to repeatedly re¬ 
synchronize the clocks of the network, the master node is often equipped with a more 
stable and precise clock than the other nodes or will have the capability to synchronize 
with an external reference time such as the one provided by the GPS. Having an accurate 
master clock is desirable in a sensor network but is not as critical as having all the clocks 
of the network synchronized with the master. 
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D. 


CLOCK MODEL 


Each device or node in a sensor network has its own clock. The time measured by 
clock A usually presents random errors, but on average it can be represented by a linear 
function of the form 

C A (t) = f A t + b A (1) 

where the term f A represents the skew and b A the offset of clock C A compared to its 
reference or master clock C B . Throughout this paper B denotes the reference node and A 
the synchronizing node. In an ideal situation with no errors or drifts, we have C(t) = t. 
The skew of a clock f A is the oscillator frequency difference between C A and C B . The 
offset b A is the time difference at a specific point in time between the two clocks. Figure 
1 is a graphical representation of a perfect reference clock combined with both a slower 
and faster slave with respective offset clocks. 



Figure 1. Clock model of a fast (C2), reference (Ci), slow clock (C3) and their 

associated equations. 

E. SOURCES OF ERROR IN CLOCK SYNCHRONIZATION 

In a perfect, ideal environment with no delayed responses, clock synchronization 
is simple. Node B sends its current time to node A, which adjusts its own clock 
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accordingly. In an actual network, many steps are required to accomplish this exchange 
of information. Associated with each of these steps is a systemic time delay which can 
be either detenninistic or stochastic. 

The following is a list of the various delays generated by the transmission of a 
timing message between two nodes [1] and illustrated using the open systems 
interconnection (OSI) model in Figure 2. The OSI model structures a communication 
system in terms of functional layers. Each layer can only communicate with the layer 
above or under it within a communication device. In order for a communication system 
to exist, two devices or nodes must share the same layered structure. 

1. Send Time 

Send time is the time delay from generating the timing message at the 
synchronizing node by its software as transfers through the application, transport and 
network layers on the sender side illustrated in Figure 2. This delay results from the 
execution of software operations and is a fixed time in most devices. 

2. Access Time 

This is the contention time to access the channel. It is highly variable and 
depends on the MAC protocol used and how busy the channel is when the node is 
attempting to access it. 

3. Transmission and Reception Time 

This is the time to send and receive the entire message by the transducer. 
This is a detenninistic delay, which depends on the size of the message sent and the 
bandwidth of the link. This is illustrated by the four green squares in Figure 2, each 
representing a packet. 
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4. Propagation Time 


The propagation is the time required for the timing message to travel 
between the two nodes. The high-latency nature of water as a communication medium 
makes detennining this delay very difficult. This problem is further explained in Chapter 
III. 


5. Encoding and Decoding Time 

This is the time required for the digital signal processor (DSP) in the 
acoustic modem to encode/decode the message to/from the received acoustic signal. This 
is a short and detenninistic delay. 

6. Receive Time 

Similar to the send time, this is the delay for the software to process the 
received timing message through the upper layers of the OSI model as illustrated in 
Figure 2. 
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Figure 2. Sources of error in estimating message latency explained using the OSI 

model. After [1]. 


To measure and eliminate some of the uncertainties associated with these delays, 
multiple messages need to be exchanged between the two nodes, adding yet more delays 
and complexity to the overall process. In some sense, clock synchronization in a sensor 
network can be regarded as the process of accounting for the effects of unknown delays 
in the overall communication process [3]. 

F. CLOCK SYNCHRONIZATION BASICS 

From Equation 1, we see that two parameters (skew and offset) need to be 
corrected to synchronize a local clock to a reference. In general, the process to execute 
this correction is achieved by transferring a set of timing messages between the two 
nodes as shown in Figure 3 and referred to as the two-way message exchange (or sender- 
receiver synchronization) [6]. 
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Figure 3. Two-way message exchange of timing messages between the synchronizing 

(A) and reference (B) nodes. After [1], 


Most authors use Figure 3 to describe the basics of clock synchronization 
protocols. The literature usually assumes a reciprocal and symmetric communication 
channel where the one-way propagation delay can be determined by averaging the total 
propagation time of the two message exchanges. Why Figure 3 can be amended to 
illustrate a time-varying propagation delay, such that the signal propagation delay 
between node A and B can be different from the one from B to A, is discussed in the next 
chapter. A solution to address this non-reciprocity is proposed in Chapter IV. 

From Figure 3, the total propagation time can be mathematically modeled as 

(T 4 -T l )-(T 3 -T 2 ) = 2p = P. (2) 

Assuming a constant propagation delay and zero or negligible skew, we can 
detennine the offset of C A from a single exchange of time-stamped messages such as 
illustrated in Figure 3. Once the propagation delays have been determined, the offset b A 
can be calculated from one of the following: 


T 2 -T\- p = b A 

( 3 ) 

T 4 -T 3 -p = —b A . 

( 4 ) 
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The offset b A from Equation 3 changes sign in Equation 4 since the offset sign is 
relative to which node it is compared with. 

By repeating the exchange of time-stamped messages, the clock synchronization 
protocol can also determine the clock skew using this fonnulation. The Time 
Synchronization for Time-Varying Propagation Underwater Network (TSVP) [7] keeps 
track of the time-stamps and calculates the clock skew using linear regression. 

Another well-known timing signaling approach is the one-way message 
dissemination method [6], The master node broadcasts its timing infonnation to all the 
nodes in its vicinity at a set time interval. A slave node collects and saves all the timing 
information in a matrix to calculate, using a least-squares estimator, the skew and offset 
of the local clock compared to the reference clock. This method assumes a fixed 
propagation delay, but the effect of this assumption is mitigated by the fact that the least- 
squares estimator in fact averages these delays. This methodology is used by the 
flooding time synchronization protocol (FTSP) [8], Although this approach requires 
multiple broadcast messages from the master node, it can synchronize a multitude of 
slave nodes simultaneously and requires minimal channel access negotiation between the 
nodes. 

How these two signaling approaches are combined and used in two different 
stages to produce a clock synchronization protocol is described in Chapter IV. 

G. LIMITATION TO CLOCK SYNCHRONIZATION 

It is possible to keep all local clocks of a network very precisely synchronized by 
simply updating them on a periodic basis. However, this scenario is unrealistic and 
would compromise the efficient use of the node’s two limited resources [6], energy and 
bandwidth, as explained in the following. 

1. Energy 

In general, devices found in sensor networks, even underwater ones, are limited in 
size and, therefore, restricted in the amount of energy they can hold. In order to 
maximize its battery and operational life, clock synchronization must be achieved while 
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preserving energy for normal operations. Since transmitting is an energy-consuming 
operation, a clock-synchronization protocol must limit the number and the length of the 
messages exchanged within the network. 

2. Bandwidth 

Underwater networks have limited bandwidth and are restricted to a relatively 
slow transmission data rate. Clock synchronization is a network support feature and not 
the main role of a network node. Limiting the bandwidth used for clock synchronization 
helps to maximize the utilization of the node for its intended main purpose. 

A balance based on the trade-offs between the desired clock synchronization 
precision and the limited network resources must ultimately be achieved to maximize the 
efficiency of the network. 
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III. VARIABILITY OF THE UNDERWATER COMMUNICATION 

CHANNEL 


The underwater acoustic channel is a challenging communication medium. It can 
be thought of as a long slim pipe, with the poor physical link quality of a mobile 
terrestrial radio channel, and the high latency of a satellite channel [9], Technologies 
developed for wireless radio communications are often used as the basis for underwater 
acoustic communications. The ocean shares with the over-the-air medium many of the 
same impairments that make wireless communications difficult such as multipath, 
scattering, noise, attenuation, and fading. The difficulties introduced by these 
characteristics are magnified in the ocean due to the high latency of the medium. 
Acoustic signals propagate five-orders-of-magnitude slower in water compared to over- 
the-air RF signals, with speeds of 1.5xl0 3 m/s and 3xl0 8 m/s, respectively. In addition, 
greater absorption at higher acoustic frequencies limits the available bandwidth, leaving 
only lower frequencies, which are susceptible to channel fading due to multipath. 

Underwater acoustic research usually takes into account the above-mentioned 
intricacies of the medium. However, most practical solutions make the assumption that 
the environment is stable over relatively short periods of time. Nevertheless, signal 
fluctuations can occur due to transceiver motion or inherent changes within the 
propagation medium. These changes are due to numerous phenomena which, depending 
on their time scale, might or might not affect the instantaneous quality of the 
communication channel. They will also have a different impact depending on whether 
the channel is in deep or shallow water. The shallow-water channel is recognized as 
being the most difficult due to its time-varying impulse response. This greatly affects the 
signal propagation delay and impacts the clock synchronization effort. 

A. TRANSDUCER MOTION 

In the underwater environment there is always some motion present in the system 
[9]. Motion can affect the impulse response’s frequency spectrum (via the Doppler 
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effect), strength of the arrival paths, or in the worst case, the complete disappearance of it 
due to one of the transceivers entering a shadow zone where no communication path 
exists between two communicating nodes. 

The Doppler effect will cause frequency shifting as well as additional frequency 
spreading. The amount of spectral shift is proportional to the ratio of the relative 
transmitter-receiver velocity to the speed of sound [9]. Since the speed of sound in the 
water is relatively slow (1500 m/s), a transducer’s motion from waves, currents, tides, 
and platform motion will be enough to impair the communication channel. Furthennore, 
unlike wireless communications systems in the air, it has been shown that assumption 
that the underwater acoustic channels are subject to wide-sense stationary uncorrelated 
scattering (WSSUS) does not always apply [10]. 

The rotation of a transducer is another dimension of motion that can affect the 
impulse response of the channel. As shown in Figure 4, in a multipath environment, the 
beam pattern of a transducer influences the directivity index for both launch angle and 
receive angle. The strongest multipath arrival may be at the start of the impulse response, 
at the end, or somewhere in between [10]. This effect is demonstrated in Figure 4, which 
illustrates a simple multipath model of an environment and its associated impulse 
response for two different transducer orientations. Rotating the transducer upward by 30 
degrees changes the relative magnitudes of the arrivals in the channel’s impulse response, 
resulting in the strongest arrival changing from the 2nd to the 1st multipath arrival. 
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Figure 4. Variation of the impulse response caused by rotating the transducer (and 
beampattern). From top left to bottom right: simulated sound-speed profile 
(isothermal), simulated multipath of the environment, beampattern when 
transducer is at 0 degrees, beampattern when transducer is at 30 degrees, impulse 
response at 0 degrees and impulse response when at 30 degrees. Note that at zero 
degrees, the second arrival has the strongest magnitude, and at 30 degrees, the 
strongest magnitude is at the first arrival. 


B. TEMPORAL VARIABILITY 

In addition to transducer motion, the impulse response will vary due to motion 
within the environment. There are various phenomena responsible for this motion, which 
fall in the two following categories [9] [11]. 
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1. Long-Term 

Changes that occur over a long time scale usually develop at a slow rate and do 
not affect the signal during its propagation through the medium. Their gradual 
occurrence means that the system can adapt to the changes in the environment. Since 
some of them can be cyclic in nature, the system can be designed to predict and adapt to 
these fluctuations. 


a. Swell 

The effect of swell is dominantly indirect, as it causes the motion of the 
transducers and affects the impulse response in many ways as explained in Section A of 
this chapter. It can also directly affect the signal reflected at the surface, inducing a 
Doppler effect on the signal frequencies based on the direction, speed and amplitude of 
the swell. 


b. Current 

Spatially, currents change the sound-speed profde of the water column. 
Temporally, the change in speed of the propagation medium intrinsically affects the 
speed of propagation of the signal. This change of speed induces a Doppler effect on the 
signal frequencies. The current creates “scintillations” (amplitude instabilities), due to 
the emergence of local turbulence in the current stream. 

c. Internal Waves 

These waves have similar effect as the current. Their inherent movement 
induces a Doppler effect on the signal frequencies. 

d. Tides 

Their period of half a day makes the change of tides to the environment 
relatively slow. They only affect communication channels located in or near a shallow- 
water area, where they modify the geometry of the multipath propagation within the 
medium. 
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e. Daily and Seasonal Temperature Variations 

Temperature variations such as the development of a warm surface layer 
affect the sound-speed profile of the water column. Again, their rate of change is very 
slow and does not affect the instantaneous characteristics of the communication channel. 

2. Short-Term 

Short-term variability is expressed in many different ways, which is the reason 
there is no consensus on a statistical characterization of the underwater acoustic channel 
[9]. Wind energy is the dominant cause of short-term instability. Its highly dynamic 
nature leaves little time for the system to adapt and can often result in the outage of the 
communication channel. 

a. Surface Wave 

The ever-moving shape of the ocean surface consistently changes the 
geometry of the surface reflection, affecting the delays of the impulse response. From 
Snell’s law, each wave crest forms an acoustic mirror with its own characteristic shape 
and focusing properties [12]. As shown in Figure 5, this change in geometry can also 
create caustics resulting in signal amplitudes much greater than those of the direct arrival 
path. The resulting impulse response can vary in time in such a manner that there is no 
clearly defined strongest arrival [9]. 



Figure 5. Qualitative description of acoustic focusing by surface wave (from [12]). 
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The movement of the wave will also, as in the case of swells, imprint a 
Doppler effect on the signal frequencies. 

b. Bubbles 

Air bubbles created and suspended at the air-water interface result in a 
significant attenuation of surface scattered signals, especially at higher signal frequencies. 
The “roughness” of the sea-surface breaking waves is directly related to the amount of 
bubbles trapped in the water [13]. 

c. Others 

There is a long list of diverse occurrences that can cause time variability 
and affect the impulse response in similar ways as explained previously. Examples of 
such occurrence are: wakes of passing ships, fish shoals, bubble screens due to rain 
showers, a fresh-water front due to river discharge, etc [10]. These examples are often 
site-specific and neglected when designing an underwater communication system. 

C. MEASUREMENT NOISE 

Measurement noise is a problem that exists in and affects practically every 
communication channel. However, interference due to noise can induce errors in the 
reception of a transmitted signal which potentially leads to errors in measuring 
propagation time-delays used in clock-synchronization protocols. Indeed, successful 
communications requires a link margin related to the signal-to-noise ratio at the receiving 
node. Like any sonar system, the variability of clock synchronization can be analyzed 
with the sonar equation. 

D. CHANNEL MODELING AND SOUNDING 

Effective single model representations of the salient propagation characteristics of 
the underwater environment have been elusive [9]. Ocean environment modeling is still 
an extremely active field of research where many aspects of its physics are still not fully 
understood. Nevertheless, the recent advances in computer processing power have made 
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it possible to better simulate the communication channel. These simulations are based on 
various models, which are not always valid for a real environment. 

In the absence of a fit-for-all underwater statistical model, channel soundings are 
a good alternative to the lack of prediction capability. They are measurements of the 
time-varying impulse response of a propagation medium [10] and replace general 
assumptions about the channel. A series of pulses is sent through the channel at a set 
interval. The receiver analyzes the received series and estimates the local channel 
characteristics from it. This procedure is an additional support cost to the system, but the 
result may greatly improve its overall efficiency. 
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IV. PROPOSED CLOCK SYNCHRONIZATION ALGORITHMS 


The following algorithms account for time variability in the underwater 
communication channel when synchronizing the clock of slave node A to reference 
master node B. As stated in the introduction, this approach is based on the TSHL 
protocol, which is explained in the following section. Our approach relaxes the 
assumption in TSHL of a stable propagation delay between the synchronizing nodes. 
Instead, we repeatedly measure the response of the channel over a period of time and 
statistically describe the communication channel to better detennine the propagation 
delay and its skew. 

Our model keeps the two-stage approach developed by TSHL, but the stages are 
implemented differently to accommodate for their respective timing message signaling 
approaches (one-way message dissemination and two-way message exchange). The 
following sections describe the original TSHL protocol in more detail followed by an 
explanation of how our approach extends each of the protocol’s stages to account for a 
communication channel characterized by a time-varying response. 

A. TSHL ACOUSTIC NETWORK PROTOCOL 

The TSHL clock synchronization protocol [1] is one of the first developed and 
published and was specifically designed for underwater or high-latency sensor networks. 
Its main feature is to split the clock synchronization process into two distinct stages based 
on the importance of correcting the skew of a clock before its offset. By first correcting 
the skew of the clock, the offset is converted from a variable function of time to a 
constant. This is necessary because otherwise the high latency of the channel makes it 
impossible to precisely correct the local clock since its offset will already have changed 
by the time the correction is applied. 

In the first stage, the skew of the clock is determined using a one-way 
dissemination method. A series of chirps (x 1 ,x 2 ,x i ,...,x n ) is sent from reference node B 
at a known and constant reference pulse interval (RPI) to the synchronizing node A. 
Linear regression based on a least-squares fit [14] is performed on the measured 
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synchronization pulse interval (SPI), which is the time interval between the arrival time 
of the pulses (y t ,y 2 ,y i ,...,y n ), relative to the known RPI. The resulting slope f A of this 

line represents the clock skew of node A. This one-way exchange of pulses is illustrated 
in Figure 6. The first stage is independent of the propagation delay which is assumed to 
be constant for the short time duration when the series of chirps is transmitted. 



Figure 6. One-way dissemination method used in the first phase of the TSHL protocol 
to correct the skew of the local clock A compared to the reference clock B. 

When stage one is completed, the skew f A between the node’s clocks has been 
detennined, and b A is the only remaining unknown value in Equation 1. The known 
skew can be corrected such that clock A’s time-line, illustrated in Figure 6, is now 
parallel to the reference clock line (horizontal). At this point, the offset b A changes from 
being a time-dependent variable to a constant since the distance between the two lines 
remains the same. 

This constant offset is corrected using a two-way exchange methodology. As 
shown in Figure 3 and assuming a constant bi-directional propagation delay, a set of 
messages is exchanged from node A to B and then from B to A. The respective node 
records and generates time-stamps, using its own clock, when the messages are 
transmitted (7, and T,) and received ( T 2 and 7’ 4 ). This timing infonnation is compiled 
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by the synchronizing node, where the total propagation time of the two messages will be 
P from Equation 2. Under the assumption that the channel is reciprocal, the one-way 
propagation time delay p k , where k is the message number, is obtained from the average 
of the two messages 

= (5) 

Knowing the one-way propagation delay from Equation 5, where p k = p , we can 
use Equations 3 or 4 to determine the clock offset b A . Both equations should yield the 
same value for the offset b A . 

Eventually, various external instabilities (listed in Chapter II) will affect the 
oscillating frequency of one or both nodes, changing the skew difference between their 
local system clocks thus requiring a repeat of the clock-synchronization protocol. 

B. STAGE 1: CHANNEL SOUNDING AND SKEW CORRECTION 

Our approach to the skew correction looks deeper into the physical attributes of 
the received chirps. When subject to a multipath environment, an acoustic modem will 
“pick” an arrival or path in accordance to a prc-dcfincd metric. In the case of the 
Teledyne Benthos acoustic modems, the most energetic path or largest peak is selected as 
the synchronizing reference datum. As explained in Chapter III, this reference arrival can 
temporally change within the impulse response of a time-varying communication 
channel. Such an environment will influence the arrivals within the impulse response 
differently, resulting in some arrivals being more stable over a period of time than others. 
The following model essentially changes the pre-defined synchronization metric used by 
the acoustic modem. The new metric is based on the structure of the impulse response 
estimated by sending a series of chirps through the channel (channel sounding). 
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Figure 7. Skew determination system diagram illustrating the transmission of the i th 

chirp of n transmitted. The response of every chirp is collected and recorded by 
the peak detector where the analysis of the ensemble begins, resulting in the 
determination of the skew of the local clock. 

The block-diagram of the proposed system is illustrated in Figure 7. 

The starting point of the process is the generation by the reference node B of a 
digital impulse train, ultimately based on its clock frequency. The length of this impulse 
train must be sufficient to provide enough chirps to the synchronization node A for the 
linear regression algorithm. Results from [1] indicate that the empirical rule of thumb for 
the optimal number of impulses required to calculate the skew is 25. The selection of the 
RPI must be based on the sampling frequency of the system and the intended 
measurement precision of the skew. The length of the whole impulse train must be such 
that there are enough samples to cover the scale of the clock skew. 

The characteristics of the pulse affect how accurately each multipath arrival can 
be distinguished from one another within the impulse response. Since the function of the 
algorithm is essentially the same as channel sounding, reference [10] identifies linear 
frequency-modulated (LFM) chirp trains and pseudorandom (PRN) binary sequences as 
the most suitable types of probe signals. A key attribute of LFM and PRN wavefonns is 
their use of a wide spectral bandwidth. Another factor contributing to the resolution of 
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the arrivals is the length of each chirp by virtue of the large time and bandwidth product 
of the probe signal. A longer chirp helps identifying the pulse within a noisy 
environment but depending on the spread of the impulse response, a long pulse can also 
introduce some superposition in the sense that successive pulses in time cannot be 
distinguished from each other. The multipath structure and distance between the two 
synchronizing nodes determines the delay spread of the impulse response and is, 
therefore, linked to establishing the length of the pulse signal. 

As the series of chirps is created, it is mixed to the carrier frequency and 
transmitted as the i th chirp ( x t (t )) of n total chirps by the reference node’s transducer. 
The environment affects the transmitted signal x i (t ), as explained in Chapter III. The 
received signal is then denoted as y t (t) in Figure 7. 

The signal y ; (f) received by the synchronizing node’s transducer is filtered, 
mixed back down to baseband, and digitized based on the sampling frequency regulated 
by the node’s local clock. 

The time delay between a signal recorded at two spatially separated sensors can 
be recovered using a correlator, even in the presence of noise [15], In this system, the 
known reference pulse x t (t) is correlated with the signal received by the synchronizing 

node y t (t) to produce R xy (r) . The distance between the peaks in R xx (t) represent the 

time delays between the multiple copies of the transmitted pulse. These originate from 
the multipath characteristic of the channel. (t) can also be used to estimate the 

channel impulse response at time t . 

The peak detector then applies a threshold to detect the location of the peaks 
present in the signal R (t). These peaks represent the signal propagation time 

difference between the different paths present in the channel at time t. 

Each arrival time is noted in Figure 7 as a { . which refers to the / h received 

multipath arrival time from the i th out of n transmitted pulses. The arrival times are then 
stored in a matrix A as illustrated in Figure 8. It is important to note that the number of 
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detected arrival can vary among the received chirp responses. Therefore, m represents 
the maximum number of arrivals detected in a single received pulse. The final size of 
matrix A is then n x m. 


j = multipath arrival times 
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Figure 8. n x m matrix A of arrival times generated by the peak detector. 

The peak detector will not necessarily detect the same number of arrivals for each 
received pulse. As explained in Chapter III, the time variability of the channel impulse 
response can affect the energy level received from a particular multipath and prevent it 
from being detected by the peak detector. Therefore, the detected arrival times present in 
a single column of Figure 8 do not necessarily represent the same multipath arrival for 
each one of the received responses. The arrivals of each multipath must, therefore, be 
properly grouped before the skew can be found. To do so, a sorting algorithm aligns each 
multipath in a single column of a matrix, keeping the received detected arrival for a 
particular pulse on the same row. 

To detennine the main multipath arrivals, a probability density function (PDF) is 
estimated using all of the detected arrival in matrix A of Figure 8. The time indexes of 
the peaks of this PDF are extracted and used to represent the main multipath, identified 
individually by the index q in Figure 9, of the channel where q = 1,2,...,<2. The 
maximum number of multipaths identified by the algorithm can be pre-set. In Figure 9, 
the s largest peaks of the PDF will be identified as the main multipath, where s = Q . 

Each main multipath is represented by a column in a new matrix B, an example of 
which is illustrated in Figure 9. The time indexes are compared with the detected arrival 
time of each pulse. If a detected arrival time corresponds to one of the determined main 
multipath arrival times, it is transferred to the new matrix on the respective row of the 
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pulse from which it was detected and the column to which multipath it is associated. If 
no arrivals are detected, a zero is inserted in the matrix as a place holder. The size of 
matrix B can be different from the size of A, depending on how much the impulse 
response of the channel varies between the transmission times of the n pulses. 


j = multipath arrival index 


i = chirp index 


IV 


a w o 


a i\ a i2 


a n\ a n2 


‘1 q 
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nq 


^Ks-n a \s 
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Figure 9. Example of matrix B generated by the PDF estimator and arrival sorting 
algorithm. B is the matrix of sorted arrival times. The number of rows in B 
remains the same as the number of rows in A (Figure 8) but the number of 
columns in B can change depending on the number of multipath arrivals identified 
by the algorithm. Zeros are inserted where no arrivals can be correlated to an 
identified multipath within a received pulse. 


Now aligned, each common ensemble of arrivals (individual columns of matrix B 
in Figure 9) is processed by a linear-regression algorithm. As mentioned in Section A of 
this chapter, the estimated time delay measured for each chirp for a common multipath 
arrival and the known RPI can be used to perfonn a linear-regression analysis. A straight 
line is fit to the data in each of the columns of B. This analysis results in two values: F j 

is the slope of the best fit line and £V represents the goodness of fit of this line, better 
defined as the linear correlation coefficient [14]. As before, the slope F- represents the 

skew between the synchronizing and reference clocks. The linear correlation coefficient 
Ej is used in the next step of the algorithm to select the most stable multipath arrival. 

The final step of the process identifies which set of arrivals is the most stable and 
reliable within the series of chirp responses analyzed. This decision is based on the 
calculated value of £V and the number of arrivals detected in a particular set. The linear 

correlation coefficient is calculated based on the standard deviations of the measured and 
expected time delays and their associated mutual covariance [14]. The value of the 
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mutual coherence can vary between 0 and 1, where 1 represents the best possible fit (i.e., 
all the points are on the line) and 0 the worst possible fit (no correlation among the 
arrivals of the ensemble). 

Time variability within a specific multipath causes the arrivals to vary in the time 
domain, resulting in a greater deviation of the arrivals from the best-of-fit line. We sort 
all of the correlation coefficients E j j = l,2,...,m and find the largest one. This most 

stable multipath arrival is chosen to be the one with the highest correlation coefficient 
(i.e., closest to 1), and this arrival has the least deviation from its best-of-fit line. 

Another metric required to confirm this most stable multipath arrival is the 
number of arrivals included in the ensemble. A low number of arrivals within a 
multipath arrival ensemble is a good indication by itself that it is not a reliable multipath 
arrival (if a column of B contains several zeros). As well, it can negatively influence the 
result of the linear regression. As an example, perfonning linear regression with only 
two arrivals produces a perfect correlation coefficient of 1. The most stable multipath is 
identified by the number of multipath detections and the associated time-variability of 
that multipath. The user must make this judgment. 

C. STAGE 2: OFFSET CORRECTION 

Our approach to correct the offset b A (as defined in Equation 1) of the skew- 
corrected clock compared to a reference is different than THSL’s as it does not assume 
that the propagation between the two nodes is the same as shown in Figure 3. The reason 
channel reciprocity is not a valid assumption in the underwater environment is explained 
in Chapters II and III. Such a situation is illustrated in Figure 10 by separately labeling 
the one-way propagation delays. 
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Figure 10. Two-way message exchange of timing messages between the synchronizing 
(A) and reference (B) nodes where the individual one-way propagation delay 
differs. In this example p A3 > P H A . After [1]. 

The definition of P in Equation 2 is modified as 

Pa,b + Pba = P • ( 6 ) 

As well. Equations 3 and 4 are modified to account for the different propagation 

delays 

T 2~ T l-PA,B= b A (V) 

P 4 ~ P 3 ~ Pb,A ~ b A ■ (^) 

The basis of our approach is the two-way message exchange used by most clock- 
synchronization protocols. Its distinction lies in its ability to differentiate the propagation 
delays for both one-way propagation channels instead of using an average based on the 
round-trip propagation time. A chirp is sent ahead of each message allowing the chirp 
response of the one-way communication channel to be detennined by the receiver. Both 
estimated impulse responses are correlated at the synchronizing node A to detennine the 
difference in propagation delay Ap between the two of them. This delay Ap is the 
difference between the time corresponding to the largest peak and zero lag. This method 
can be used because the modem always synchronizes on the largest amplitude arrival. 
The offset of the clock can be calculated using this propagation delay difference and the 
four time-stamps generated by the two-way message exchange. 
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An overall system diagram of this adapted two-way message exchange is 
presented in Figure 11 and is further explained below. 



Figure 11. Adapted two-way message exchange for calculating the offset of a skew- 

corrected local clock compared to a reference. 

The synchronizing node A starts the offset correction process by sending a new 
pulse ahead of the offset calculation request message to reference node B. The time of 
the start of the transmission is time-stamped T x and recorded by node A. 

As it propagates through the channel, the pulse is transformed as explained in 
Chapter III. The yet-to-be-detennined propagation delay is denoted as p A3 in Figure 11. 
Time-stamp T 2 is recorded by the receiver upon first detecting the arrival of the pulse 
sent by node A. From the pulse response, the instantaneous one-way impulse response 
from node A to B is determined using a matched filter and identifying the peaks of the 
resulting time-series. This impulse response is required later in the algorithm at the 
synchronizing node. To reduce the bandwidth required to send it across the channel for 
processing, the estimated impulse response is compressed by retaining its main 
characteristics, which are the amplitudes of the largest peaks and their associated time 
indices. These characteristics are represented as H A B in Figure 11. 
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The reference node then sends the same new reference pulse before transmitting 
the second message of the exchange. Time-stamp T 3 is recorded concurrently to the 
beginning of the transmission of the second message. It carries the two time-stamps 
generated by node B, T 2 and T 3 , as well as the compressed impulse response information 
H AB recorded from the probe signal preceding the first message. 

The second message propagates through the channel from node B to A. The 
message is transfonned by the channel as discussed in Chapter III. The yet to be 
detennined propagation delay is denoted as p BA in Figure 11. Time-stamp T 4 is 

generated by the receiver upon detecting the arrival of the chirp sent by node B. From 
this received chirp, the instantaneous one-way impulse response from node B to A 
(H BA ) is estimated and its infonnation compressed using the same process explained in 
step 3 of this algorithm. From their respective recorded characteristics, the main arrivals 
of both estimated impulse responses are reconstructed in the time-domain from H A3 and 
H B A . Their cross-correlation allows for the time-delay estimation between the two 

signals. Based on the assumption that the modem will synchronize with the strongest 
arrival, the difference in propagation delay A p (Equation 9) will be the time index 
difference between the highest peak of the cross-correlation and its zero lag: 

Pa,b-Pb,a = A P = z *- ( 9 ) 

The total propagation delay P of both messages and the time-delay difference 
A p can be combined to evaluate one of the two propagation delays ( p BA or p AB ). 
Adding Equations 6 and 9 together, we get 

P + Ap = 2p APi . (10) 

The one-way propagation delay p AB can be determined from Equation 10: 


P + Ap 

Z — Pab ■ 


( 11 ) 


The second propagation delay (in this case p BA ) can be evaluated from 
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P Pab Pb,a • 


( 12 ) 


Either one of the one-way propagation delay can be used in Equation 11 as A p 
can be positive or negative. Whichever one-way propagation is evaluated in Equation 11, 
the remaining one can be detennined using Equation 12. 

The last step of the algorithm consists of evaluating the offset b A between the two 
nodes’ clock using Equations 7 or 8. The result of both equations should be the same. 

D. TIME REVERSAL APPROACH 

A time reversal approach was investigated as a way to evaluate the individual 
one-way propagation delays of the two-way message exchanged methodology. This 
approach initially appeared very promising since one of its benefits is the elimination of 
the multipath effect of the communication medium. However, one of the assumptions of 
time-reversal is the requirement of reciprocity within the channel [16]. The basis of our 
algorithm is the existence of different one-way propagation delays resulting in a non¬ 
reciprocal channel. Therefore, a time reversal approach was rejected for our problem. 
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V. COMPUTER SIMULATIONS STUDIES 


We implemented the approach described in Chapter IV in Matlab. The Bellhop 
and VirTEX underwater acoustic propagation models (both compatible with Matlab) are 
used to simulate a simple underwater acoustic communication channel. The aim of the 
simulation is to prove our conceptual approach prior to conducting an in-water 
experiment. The environmental models, the simulated environment’s parameters, the 
step-wise results, and the clock-correction stages (skew and offset) are explained in the 
following sections. 

A. BELLHOP 

Bellhop [17] is a ray-tracing program developed by Michael B. Porter and 
distributed as part of the Acoustic Toolbox software package. Based on ray theory, it can 
produce two-dimensional acoustic ray tracing for a given channel geometry and 
underwater environment. Various options can be selected and inputs given to the model 
to better define the environment and provide more realistic simulation. Once the channel 
geometry and environment are defined, Bellhop provides various graphical and numerical 
outputs. 

B. VIRTEX 

The VirTEX family of algorithms was designed by Peterson and Porter [18]. The 
algorithms build on the results obtained by a Bellhop simulation and add the effects of 
multipath and Doppler introduced by environmental motion. The algorithm computes the 
time series that would be observed at a hypothetical receiver location. 

C. DESCRIPTION OF THE SIMULATED ENVIRONMENT 

Our simulation environment is very simple. The aim of our simulation is not to 
model a specific environment but to test the theory of our algorithm. The environment, 
illustrated in Figure 12, consists of a transmitter and a receiver, both situated at a depth of 
100 m and separated by a distance of 1 km. The sound-speed profile is isothermal for 
simplicity. Time variability is introduced in the environment by randomly changing two 
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parameters: the orientation of the transmitter and the phase of a simulated wave on the 
sea surface. As explained in Chapter III, the orientation of the beam pattern can change 
the path which has the most energy. As well, surface waves can alter the delay-spread 
characteristics of the impulse response. The surface wave is modeled as a sinusoidal 
wave with phase randomly changed throughout the simulation to provide a non-stationary 
surface boundary. 


Normally, both nodes of the system will operate with the same sampling 
frequency. To simulate a skew difference between the two clocks, the sampling 
frequency of the synchronizing node is slightly altered when recording the incoming 
chirp time-series. This artificial method creates a kn own skew between the two nodes’ 
clocks. For this simulation, the sampling frequency at the receiver is reduced from 
48,000 samples/sec to 47,885 samples/sec. This reduction in sampling frequency results 
in a skew of 417 ppm. 
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Figure 12. Characteristics of the simulated environment. From the top left corner: the 
sound-speed profile; the geometry of the environment and its associated ray paths; 
the impulse response of the communication channel. 
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D. STEP-BY-STEP SIMULATION RESULTS 

This section presents the formal results of the simulations. Figures illustrating the 
different stages of the algorithm are also presented to help visualize some of the 
explanations from Chapter IV. 

1. Phase 1: Skew Correction 

A 5-ms LFM chirp based on a 5-kHz bandwidth is created and designated as the 
reference chirp. The chirp is sent through the channel 20 times at a RPI of one per 
second. The orientation of the transmitter and phase of the surface wave is randomly 
changed prior to each transmission and remains the same during the transmission of a 
single chirp. The simulation is run with and without a skew difference. 

The receiver, knowing the sampling frequency, RPI and number of pulses sent, 
separates the pulses into a matrix of 20 chirps as illustrated in Figure 13. 
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Chirp responses 



Figure 13. Ensemble of 20 chirp responses with no clock skew. 


Without any signal processing, time variability is already noticeable in Figure 13. 
Let n denote the number of probing pulses (chirps) used in any given experiment, n = 20 
is then the number of measurement waveforms displayed in Figure 13. We denote these 

n waveforms by { V,(0}" =l • 

We compute the cross correlation R (t) between the reference chirp x(t ) and 
each of the measurement signals y,(t) to create the ensemble of cross correlations 

■[/?„, (f)j i plotted in Figure 14. 
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Cross correlated responses from chirps 



Figure 14. Ensemble of 20 cross-correlated received chirp responses with no clock skew 

present. 

Note that the signals in Figures 13 and 14 contain m = 5 multipath arrivals that 
make up the chirp response. Let j denote the “arrival index,” or the index on the m 

r -I m 

multipath arrivals that occur in a given single measurement 

waveform, j = 1,2,..., m . 

The peaks displayed in Figure 15 are estimated from the cross-correlation signals 
displayed in Figure 14. First, the individual arrivals in the cross-correlation waveforms 
are selected to be the multipath arrivals that exceed a small threshold chosen by the user. 
For each of the arrivals in a given correlation wavefonn we find the time x" 

corresponding to the time at which a peak occurs by choosing the time lag t such that 

^,.(0 = max {^( T )}. 
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Figure 15. The peak times for an experiment in which clock skew is present: The peak 
times Uji derived from the ensemble of n = 20 measurement signals depicted in 

Figure 14 are plotted in this figure. Note that each of the m - 5 arrival wavelets 

cij, j = 1,2,..., m has a corresponding ensemble U } = [n ;i J , j = 1,2,..., m 

associated with it. 

We see in Figure 14 the ensemble of n waveforms and in Figure 15 the ensemble 
of peaks of waveforms that result from the n = 20 probing pulses from the transmitter. 
After peak detection, we see in Figures 15 and 16 that each multipath arrival has a 
corresponding ensemble of n = 20 peaks associated with it. The peak arrival times shift 
over time due to the skew between the transmitter and receiver clocks. An example of an 
ensemble of waveforms with zero skew is depicted in Figure 16. 
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Figure 16. The peak times for an experiment in which zero clock skew is present: The 
peak times n ;7 derived from the ensemble of n = 20 measurement signals depicted 

in Figure 14 are plotted in this figure. Note that each of the m = 5 arrival 
wavelets a., j = 1,2,..., m has a corresponding ensemble 


Uj = {«, | , j = l,2,...,m associated with it. 


Let i denote the peak index corresponding to each of the n peaks that occur for a 
single multipath arrival in a given ensemble constructed from the data displayed in 
Figures 15 and 16. Because n probe pulses create an ensemble of n waveforms, there 
are n peaks in an ensemble corresponding to the j th multipath arrival in Figures 15 and 

16, i = l,2,...,n. We denote the z th peak time for the y' th multipath arrival by u j{ . We 

denote the ensemble of peak times for the/ h multipath arrival by U ; = : , 

j = 1,2,..., m . 

For each multipath arrival a ., j = 1,2,..., m , we estimate a PDF of peak arrival 
times using the ensemble Uj = {n ; , j , j = 1,2,..., m . We use a kernel density estimator 
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with a Gaussian kernel to build the PDF estimates [19]. Once all m PDFs are estimated 
individually, we plot them on a single graph as depicted in Figure 17. We see that in 
Figures 15, 16 and 17 there are m — 5 multipath arrivals. 

Individual PDFs of the arrivals for each path (blue) and impulse response of the channel (red) 



Figure 17. The scaled simulated impulse response (a set of Kronecker delta functions) of 
the channel is plotted in red. The estimated PDFs of the individual arrival 
wavelets are plotted on a single graph and appear in blue. Note that the delta 
functions occur at the times at which the chirp pulses begin, so the centers of the 
PDF estimates are shifted slightly to the right of delta functions. 

The shape of each PDF, as illustrated in Figure 17, is representative of the 
stability of each signal propagation path. The direct paths and bottom-bounce paths are 
not affected by the time variability of the environment, which is only present at the water 
surface interface. These unaffected paths are characterized by sharp and narrow PDFs 
whereas the surface-bounce paths are wide and low. To demonstrate the time-variability 
of the environment, the non-Doppler affected estimated impulse response, as shown in 
Figure 12, is superimposed on Figure 17. It is important to note that the amplitudes of 
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the PDFs are not related to the amount of energy propagated through the associated path. 
The third arrival shown in Figure 17 is a good example of the case in which the weakest 
signal path is very stable. 


Individual PDFs of the arrivals for each path (blue) and impulse response of the channel (red) 
x io -3 zero clock skew present 



Figure 18. The scaled simulated impulse response (a set of Kronecker delta functions) of 
the channel is plotted in red. The estimated PDFs of the individual arrival 
wavelets are plotted on a single graph and appear in blue. Note that the delta 
functions now occur in the middle of the chirp pulses, as zero skew is present. 

Figure 18 is a replica of Figure 17 but with zero skew present. Since the 1st and 
3rd arrivals are not affected by the time variability of the channel and zero skew is 
present, their respective PDFs are now represented by Kronecker delta functions. In 
contrast to Figure 17, the delta functions now occur in the middle of the chirp pulses. 

A least-squares linear regression algorithm is applied to each ensemble of arrivals. 
The result of this linear regression is the clock skew difference F J , from which the 

coefficient of determination E, can be calculated as explained in Chapter IV. Since the 

tolerance of a clock is usually measured in ppm, Equation 13 can be used to convert the 
skew to Fpj to compare the measured value with the published tolerance of a clock: 
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\,QQ0,000(F j -\)=Fp j . 


( 13 ) 


The simulation results for the first phase of the clock-synchronization algorithm 
are illustrated in Table 1. The algorithm chooses the first path as the most reliable 
(shaded arrival in Table 1) as it has the highest correlation coefficient (perfect in this 
case). This is expected, since both the first and third arrivals do not interact with the sea 
surface and, therefore, no time variability is introduced in their associated signal 
propagation path (direct and single bottom bounce) as illustrated in Figure 12. 


Table 1. Skew correction simulation results. The shaded result corresponds to the most 
stable arrival set as selected by the algorithm (the set with the largest correlation 

coefficient). 


Simulation 
and set clock 
skew 

Arrival 

Ensemble 

Wj ) 

Skew ( Fj) 

Correlation coefficient 
(■ E J ) 

Skew in 
ppm 

( F Pj ) 

#1 

No clock 
skew 
(0 ppm) 

1 

1.000000000000000 

1.000000000000000 

0 

2 

1.000087954260652 

0.999999926902282 

88 

3 

1.000000000000000 

1.000000000000000 

0 

4 

1.000289301378446 

0.999999011969708 

289 

5 

1.000683255012531 

0.999995976887908 

683 

#2 

With clock 
skew 

(417 ppm) 

1 

1.000416666666667 

1.000000000000000 

417 

2 

1.000504620927318 

0.999999926963154 

504 

3 

1.000416666666667 

1.000000000000000 

417 

4 

1.000705968045113 

0.999999012792314 

705 

5 

1.001099921679198 

0.999995980236109 

1099 


After artificially introducing a clock skew in the simulation, the algorithm 


successfully detennines the skew as 417 ppm, which is the value expected from the 
reduction in sampling frequency at the synchronizing node (48,000 samples/sec to 47,885 
samples/sec). 

This simulation demonstrates that time variability can affect the measurement of 
the skew between two clocks. The skew values calculated from the ensembles of the 
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multipath arrival times affected by time variability varies greatly. In the worst case 
(ensemble #5) the skew value in ppm is nearly three times the real skew value. 

2. Stage 2: Offset Correction 

The same environment and reference chirp as in stage 1 are used to simulate our 
approach to the clock-offset correction. The simulation is designed to follow the steps 
explained in Chapter IV and illustrated in Figure 11. It focuses on the calculation of the 
difference Ap in propagation time between the two messages sent in the two-way timing 
message scheme. Once this parameter is determined, the other unknown values 
(P’Pab’Pba an d b A ) can be found using the four recorded time-stamps and Equations 2, 

11 and 12. An example of how these parameters can be calculated with these equations 
is illustrated in Figure 19. 

Synchronizing Reference 

Node A Node B 

Ti= 50 

r 4 = 73 

i Ap= ~ 2 j 

Figure 19. Example of time-stamp exchange to detennine the clock offset b A between 

nodes A and B. 

From Equation 2, the total propagation time P is calculated as 

(T 4 -T 1 )-(T 3 -T 2 ) = (63- 50)-(40 - 35) = 18 = P. (14) 

Then, using Equation 11, and knowing P and Ap , we get 




46 



= 10 = P a ,b- 


(15) 


P - Ap 
2 


18-(-2) 


The last propagation delay p B A can be calculated from 

P ~ Pas = 1 8 —10 = 8 = p BA . (16) 

Finally the offset b A can be detennined using either Equations 7 or 8. 

To test our approach, the simulation must create an environment where a Ap will 
exist between the two message propagation delays. Assuming that the modem (node) 
picks the strongest arrival to synchronize with, we see that the simulation must vary the 
location of this arrival within the impulse response of the two message transmissions to 
create a Ap . The parameters used in Figure 4, where the transmitter was tilted from 0 to 
30 degrees, were used again for this simulation since they illustrate the desired temporal 
displacement of the strongest arrival. The steps of the simulation, illustrating the various 
signal processing steps, are described in this section. 

The synchronizing node A initializes the clock offset correction phase following 
the determination of the clock skew in stage 1. It sends a reference chirp and an offset 
correction request message to the reference node B while also creating and storing a time 
stamp 7, when the transmission is initiated. 

The received chirp response, illustrated in Figure 20, is processed, where its main 
characteristics (peak amplitudes and corresponding time index) are extracted in order to 
limit the amount of information transmitted across the channel. Time stamps T 2 and T 3 
are respectively created by the reference node when it receives the first message and 
transmits the second one. The difference between these two time stamps represents the 
processing or handling of the request by the reference node. The reference node 
compiles the two time-stamps and the characteristics of the analyzed estimated impulse 
response into a message, which is sent back to the synchronizing node following the 
transmission of another chirp. 
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Figure 20. LFM chirp response received by the reference node B. 

The synchronizing node receives the chirp response and processes it in the same 
manner as did the reference node. This results in a basic estimated impulse response for 
each one-way transmission between the two nodes, illustrated in Figure 21. 
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Figure 21. 


Simplified estimated impulse responses from the two transmissions 


These two estimated impulse responses are cross-correlated resulting in the cross 
correlation signal plotted in Figure 22. 
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Figure 22. Cross-correlation of both simplified impulse responses. The red line indicates 
the middle of the cross-correlation. The distance between the highest peak and the 
middle of the cross-correlation represents A p . 

The delay between the highest peak and the center of the cross-correlation (time 
index 0) is the difference in propagation delay A p between the two messages sent. The 
measured time series simulated by VirTEX is referenced in the time domain to the 
beginning of the transmission, as shown in Figure 20, which allows the determination of 
the “true” propagation delays p B A , p AB and A p using Equation 9. 
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Table 2. True and calculated values of the difference in propagation delays between the 

two messages during the two-way exchange. 


True P AS 

0.6664 sec 

True P b ,a 

0.6536 sec 

True Ap 

-0.0128 sec 

Measured Ap from cross-correlation 

-0.0128 sec 


As stated at the beginning of this section, once Ap is detennined, the clock offset 
is easily calculated. The simulation proves the theory of our approach since it 
successfully and precisely determines the difference in propagation delays and clock 


offset. 
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VI. DEL MONTE LAKE EXPERIMENT 


A. OBJECTIVE 

The aim of the experiment is to test and prove the algorithms described in Chapter 
IV in an actual, time-variant, underwater environment. Only the first stage (skew 
correction) is tested during the experiment. The second stage requires resources not 
available to this project, such as access to the firmware of the Teledyne Benthos acoustic 
modems. 

Del Monte Lake, situated on the campus of the Naval Postgraduate School is 
selected for its accessibility and minimal environmental studies requirements and 
approvals. The general concept of the experiment consists of transmitting and recording 
a series of chirps between two acoustic modems. The recorded data are downloaded 
from an on-board digital recorder and processed separately. During the post-experiment 
processing, the arrivals of chirp responses are tracked and processed to determine the 
skew difference between the clocks of the two acoustic modems (nodes). The lake 
environment on the day of the experiment, the equipment used, the test procedures, 
selection of parameters and an evaluation of the results are explained in more detail in the 
following section. 

B. ENVIRONMENTAL DESCRIPTION 

The experiment was conducted during the afternoon of 2 May 2012 on the Del 
Monte lake. A satellite picture [20] of the area is illustrated in Figure 23. 
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Figure 23. Satellite picture of Del Monte Lake and locations of the acoustic modems 

during the experiment (from [20]). 


It was a warm sunny day with reported winds of 24 km/h [21]. The wind was 
strong enough to generate waves on the surface of the lake and made it difficult to 
maintain a stationary location with the small aluminium boat. From a previous 
experiment on the lake, the depth of the water was found to be approximately 2 m over 
most of the lake. The bottom of the lake consists of thick mud. The water temperature 
was not measured but was assumed to be isothermal through the water column due to the 
very shallow depth. The water fountain (not shown on Figure 23) was turned off to 
reduce the amount of noise present in the lake. 


C. EQUIPMENT DESCRIPTION 

1. Aluminum Boat and Motor 

A 13 foot aluminum boat equipped with a Minn Kota ENDURA Pro C2 transom- 
mounted trolling freshwater electric motor was used to conduct the experiment. It was 
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used to deploy the acoustic modem and served as a working platform on which the deck 
box and laptop were operated to perform the experiment. 

2. Teledyne Benthos SM-75 SMART Modem 

A Teledyne Benthos SM-75 Smart Modem Acoustic Release Technology 
(SMART) [5] was refurbished and fitted with an onboard digital recorder, as illustrated in 
Figure 24. 



Figure 24. Teledyne Benthos SM-75 SMART modem (from [5]). 

The modem operates in the low-frequency (LF) range of 9 -14 kHz and uses a 
sampling frequency of 10240 samples/sec. The digital recorder is equipped with a 4-GB 
Secure Digital High Capacity (SDHC) card, which records the baseband received signal 
by the modem in minute-long .wav or .pern files. The gain of the received signal is 
continually adjusted, using an automatic gain controller (AGC), to keep the signal level 
constant in the recording. The AGC information is also recorded as part of the baseband 
signal by the digital recorder. A command is sent to the SM-75 to start and stop the 
digital recording, which allows the recording of selected desired events. Once recovered, 
the SM-75 can be dismantled and the SDHC card containing the recorded data can be 
removed. These recordings are then uploaded to a stand-alone computer for post- 
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processing and analysis. For better results, Teledyne Benthos provides a Matlab routine 
which uses the AGC infonnation of a recording to restore the signal to its original level 
and continuity. 

The LF beam pattern of the SM-75 is illustrated in Figure 25. This beam pattern 
is used in Chapter II to demonstrate the effect of a transducer rotation on the impulse 
response of the channel. 
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- 150.0 
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Figure 25. SM-75 LF (9-14 kHz) beam pattern (from [5]). 

During the experiment the SM-75 was deployed and anchored at station A, as 
shown in Figure 23, at a depth of approximately 1 m. It functioned as the synchronizing 
node B, recording the series of chirp response transmitted by the reference node A (deck 
box). 

3. UDB-9400 Universal Deck Box 

As illustrated in Figure 26, the deck box is an acoustic modem where the 

transducer and modem electronics have been separated. Several interfaces, including a 
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touch screen display, are combined with the modem in a portable case, which allows a 
user to perform network administration within the network. The transducer is linked to 
the portable case via a long cable and can be deployed over the side of a boat or platform. 



Figure 26. UDB-9400 Universal Deck Box, including its deployable transducer (left) 

(from [5]). 

The modem part of the deck box uses the same communication parameters used 
by the SM-75 as described in the previous section. The deck box is not equipped with a 
digital recorder but still has a limited flash memory module, where files can be uploaded. 
For this experiment, .wav audio files of the various series of chirps were uploaded to the 
memory from which they could be played. Instead of using the touch screen to input 
modem commands, the deck box was connected via serial port (RS-232) to a laptop 
running the PROCOMM software. PROCOMM provides an easier user interface and 
allows the creation of script files that made the conduct of the experiment easier and 
faster to complete. The same laptop can also be interfaced to the SM-75 or acoustic 
modem. 

The deck box was used as the reference node A during the experiment. Various 
series of chirps were transmitted at a distance of 2 m (station B on Figure 23) and 50 m 
(station C on Figure 23). The transducer was attached half-way down a 3 m long metal 
pole which was lowered and embedded in the mud at the bottom of the lake. The 
transducer was kept stable this way in the middle of the water depth at approximately 1 
m. 
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D. EXPERIMENTAL PARAMETERS 


The parameters of the experiment had to be chosen to accommodate the 
characteristics of the test environment. Underwater acoustic modems are not normally 
deployed in a shallow fresh water lake and limited information was available about any 
previous experiments in such environments. To mitigate the risk associated with this 
unknown environment, various parameters were tested during the experiment. These 
parameters are explained in the next sections. 

1. Distance Between the Nodes 

The shallow depth of the lake makes the communication channel prone to 
multipath. The simultaneous arrival of many ray paths at the same time can make the 
retrieval of any distinct signal impossible due to temporal superposition. For a 2 m water 
depth, the maximum propagation delay between the direct and surface or bottom bounce 
paths is obtained by separating the nodes by 2 m. Assuming an isothermal water column 
and a constant sound of speed of 1530 m/s, the modeled time-delay difference between 
the two main paths is calculated to be 0.0005 s compared to 0.0001 s when the nodes are 
separated by 50 m. These two distances were used during the experiment. The ray paths 
and impulse response for nodes separated, respectively, by 2 m and 50 m were modeled 
using Bellhop and are shown in Figure 27. 
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Figure 27. Bellhop model of the ray paths and impulse responses in Del Monte Lake with 

2 m and 50 m node separation. 


2. Probe Signal Characteristics 

The design parameters of the probe signals are important as they determine how 
well each ray path can be distinguished at the receiver. It is generally accepted that the 
LFM chirp is a good choice as a probe signal [10]. The experiment used both LFM 
chirps and CW pulses in order to give some measure of the advantage of using the LFM 
chirp. The frequency of the chirp swept over the range of 0 to 5000 Hz which was later 
modulated by the modem to its operational bandwidth of 9 to 14 kHz. 

The length of the probe signal affects the amount of temporal superposition at the 
receiver. As an example, when the nodes are separated by 2 m, the time-delay difference 
between the two main paths is such that a 0.5 ms probe signal will not be superimposed. 
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At a sampling frequency of 10240 Hz, a 0.5 ms probe signal corresponds to only 6 
samples in length. Based on the fact that the ambiguity function of the chirp is inversely 
proportional to its length [22], a longer chirp should provide better time resolution. 
During the experiment, chirps of 0.5, 1, 2, 5, 10, 25 and 50 ms duration were used and 
compared in the post-processing analysis. 

The pulse repetition interval was selected and remained at 250 ms throughout the 
experiment. This value was arbitrarily selected as it seemed short enough to respect the 
requirement to minimize the utilization of the network bandwidth when perfonning a 
clock synchronization procedure (as explained in Chapter II). Another factor that was 
considered was that the interval be short enough so that the AGC does not reset in- 
between each probe transmission. 

The final parameter was the number of pulses sent in a given series. This was 
selected to be 30 based on the experimental result of Syed [1]. During that experiment, it 
was detennined that 25 pulses was the optimum number of pulses before observing a 
diminishing benefit in increasing the number of pulses. Since Syed’s experiment was 
conducted over the air, our number of pulses was increased to 30 to ensure that 25 of 
them would be received in the less favorable underwater environment. 

E. TEST PROCEDURES 

The experiment was executed in three stages, which translate into the preparation, 
the acquisition of measurements, and the analysis. The procedures associated with each 
stage are described in this section. 

1. Preparation 

In preparation for the experiment the above-stated experimental parameters were 
chosen. A series of 30 identical pulses was generated using Matlab and then converted to 
a .wav file format. Each file was uploaded to the deck box’s internal flash memory via 
the serial port using a standard laptop computer. Once loaded in the deck box’s memory, 
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the .wav file could be played through the modem and its transducer. A PROCOMM 
script file was written to execute the experiment more efficiently and with fewer manual 
errors. 

Before proceeding to the lake, an over-the-air experiment was conducted in the 
Seaweb Lab to confirm that the equipment worked as planned. The original script called 
for the digital recorder to be turned on and off between the transmission of each series of 
pulses, which sometimes caused the SM-75 to reboot randomly. The script was modified 
so that the digital recorder was turned on at the beginning of the experiment and only 
turned off once every series of pulses was transmitted. This solution was effective and 
only one measurement was lost due to a SM-75 reboot. 

2. Experimental Measurements 

The equipment was loaded into the aluminum boat and transported to the 
deployment stations as shown in Figure 23. The SM-75 was deployed and its location 
recorded using a handheld GPS. The 2-m separation between the SM-75 and the deck 
box transducer was confirmed by the handheld GPS and by the ranging routine included 
in the onboard software of the Teledyne Benthos modems. Every pre-recorded series of 
probe signals was transmitted twice at both the 2-m (experiment #1 and #2) and 50-m 
(experiment #3 and #4) stations, for a total of four experiments. No additional 
experiments were perfonned due to the drained battery of the laptop which was essential 
to the conduct of the experiment. The SM-75 was recovered and the equipment was 
returned to the Seaweb Lab. The recordings from the SM-75 were transferred to a laptop 
for analysis. 


3. Post-experiment Analysis 

A log file was created along with each recording by the digital recorder. The log 
files identify every .wav file of duration of one-minute long or less associated with it. 
Time-stamps were also generated by the modem’s internal clock and annotated in the log 
marking the beginning and ending of each files. As stated earlier in this chapter, 
Teledyne Benthos provides a Matlab routine that reconstructs the recording based on the 
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log file. With this routine, each recording was processed and the original non-AGC time 
series were then available in a Matlab format. 

Each recording contained multiple transmissions of pulse series. These series 
were easily identifiable, as shown in Figure 28, and were individually extracted to be 
processed separately. 



Figure 28. AGC-corrected recording where all series of probe signals were received 

when the nodes were separated by 50 m (experiment #3). 

Each individual series of probe signals was processed separately using the same 
algorithms previously described in Chapter IV and V. 
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F. RESULTS 

The cumulative results of the analysis are illustrated in Table 3. Following Table 
3, a selection of these results is expanded to demonstrate some of the noteworthy 
findings. 


Table 3. Estimated skew (in ppm) and its related error for each different series of probe 
signals transmitted during each experiment. An average of the error (in ppm) for 
each different series of probe signals is also included. 


Estimated skew and error in ppm for each experiment (all values are in ppm) 

Probe Signal 

Experiment 
#1 (2m) 

Experiment 
#2 (2m) 

Experiment 
#3 (50m) 

Experiment 
#4 (50m) 

Average 

Error 

0.5ms LFM 

-42 ± 10 

-58 ± 15 

-33 ± 14 

-28 ±6 

±11.25 

1ms LFM 

-42 ± 10 

-21 ±25 

4 ± 16 

-55 ± 14 

±16.25 

2ms LFM 

-28 ± 12 

-53 ± 17 

-45 ±3 

-35 ±4 

±9 

5msLFM 

-38 ± 11 

-34 ±11 

-33 ±11 

-29 ±4 

±9.25 

lOmsLFM 

- 

-44 ± 50 

-20 ± 14 

-28 ±5 

±23 

25msLFM 

-41 ± 14 

-54 ± 22 

-31 ±6 

-51 ± 10 

±13 

50ms LFM 

60 ±27 

-32 ± 27 

-34 ±11 

11 ± 15 

±20 

0.5ms CW 

11 ±34 

5 ±46 

-56 ±25 

-48 ±3 

±27 

1ms CW 

-46 ±8 

5 ±33 

-8 ±31 

-15 ± 15 

±21.75 

2ms CW 

70 ±38 

-44 ±4 

-36 ± 20 

-40 ± 27 

±22.25 

5ms CW 

-42 ± 33 

9 ±35 

-13 ± 16 

-59 ±21 

±26.25 

10ms CW 

-30 ±32 

-39 ±39 

-28 ±4 

-16 ±26 

±25.25 


The digital recorder rebooted during the playing of the 10ms LFM series of pulses 


in experiment #1, therefore, no results are shown for this event in Table 3. 

An initial glance seems to demonstrate a high variance among the results. To 
better analyze this variance, the uncertainty of each skew measurement is calculated [14]. 
The following Figures and Tables examine the source(s) of this variance and provide a 
higher degree of confidence in the experimental results. 

From Table 3, the 2 ms LFM chirp is found to be the superior probe for this 
environment based on its lower average error of ±9 ppm. In Figures 29 and 30, 
respectively, the first 15 matched-filtered responses of the 2ms chirp when tested at 2 m 
(experiment #2) and 50 m (experiment #3) are shown. These figures are used as a 
comparison with other chirps that were identified to be less effective. 
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Figure 29. 


Figure 30. 


First 15, 2ms chirps (experiment #2) 



Time sample index 


1200 15 


Chirp index (i) 


First fifteen 2 ms long LFM chirp responses at 2 m (experiment #2). 


First 15, 2ms LFM chirps (experiment #3) 
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Looking at Figures 29 and 30, we see that the first detected and most likely the 
direct path is easily identifiable and looks fairly stable over the measurements. The 
following arrivals are not as easy to distinguish but a few peaks are estimated to be 
associated with a multipath arrival by the algorithm as shown in Tables 4 and 5. The data 
presented in these tables are the time index where a peak is detected. Perfonning linear 
regression on each of these paths allows the determination of the clock skew and its 
correlation coefficient, which is used to identify the most stable path. The shaded column 
highlights the most stable path based on its correlation coefficient. In both cases, the first 
path has the highest correlation coefficient, which is consistent with Figures 29 and 30. 
The first detected path clearly shows a slow shift of the detected peak, indicating the 
presence of skew between the clocks in the two nodes. Respectively following Tables 4 
and 5 are Figures 31 and 32, which illustrate the same data presented in the tables. These 
figures clearly show the time variability of the multipath arrivals. 


Table 4. Detailed results of the 2 ms LFM pulse received during experiment #2, including 
the correlation coefficient, clock skew (in absolute and ppm values) and its 

associated error. 


Analysis of 2 ms LFM pulse received during experiment #2 (2 m separation) 

Pulse Index 

Time index of detected peaks for each multipath 

Detected 
path #1 

Detected 
path #2 

Detected 
path #3 

Detected 
path #4 

Detected 
path #5 

1 

362 

389 

530 

557 

599 

2 

363 

394 

535 

0 

0 

3 

363 

389 

540 

566 

600 

4 

363 

395 

528 

0 

0 

5 

363 

0 

531 

0 

611 

6 

363 

395 

531 

564 

599 

7 

363 

390 

535 

565 

599 

8 

364 

395 

531 

557 

599 

9 

364 

390 

532 

564 

602 

10 

364 

395 

532 

564 

603 

11 

364 

395 

532 

559 

0 

12 

354 

381 

532 

563 

603 

13 

364 

395 

532 

0 

0 
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Pulse Index 

Time index of detected peaks for each multipath 

Detected 
path #1 

Detected 
path #2 

Detected 
path #3 

Detected 
path #4 

Detected 
path #5 

14 

361 

393 

532 

562 

589 

15 

364 

396 

532 

566 

593 

16 

364 

390 

532 

0 

0 

17 

365 

396 

532 

566 

600 

18 

364 

395 

533 

0 

609 

19 

365 

391 

532 

566 

600 

20 

365 

391 

532 

559 

604 

21 

365 

396 

0 

0 

0 

22 

365 

395 

533 

0 

601 

23 

365 

396 

533 

0 

604 

24 

365 

391 

533 

564 

595 

25 

365 

396 

534 

0 

609 

26 

366 

396 

531 

565 

0 

27 

366 

396 

0 

0 

604 

28 

366 

394 

0 

0 

0 

29 

366 

392 

534 

0 

0 

Correlation 
Coefficient ( E f ) 

0.99999999 

0.99999997 

0.99999998 

0.99999996 

0.99999992 

Skew (Fj ) 

1.00005291 

1.00004439 

1.00000411 

1.00004903 

1.00003269 

Skew (Fpj ) 

-53 

-44 

-4 

-49 

-33 

Skew error 

±17 

±30 

±21 

±44 

±61 
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Detected Peaks of 2 ms LFM chirp from experiment #2 (2m) 
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Figure 31. The peak times for an experiment in which clock skew is present: The peak 
times u jt derived from the ensemble of n = 30 measurement signals depicted in 

Figure 29 are plotted in this figure. Note that each of the m = 6 multipath arrivals 

a.j, j = 1,2 ,..m has a corresponding ensemble U j = j ,j = l,2,...,m 

associated with it. 


Table 5. Detailed results of the 2 ms LFM pulse received during experiment #3, including 
the correlation coefficient, clock skew (in absolute and ppm values) and its 

associated error. 


Analysis of 2 ms LFM pulse received during experiment #3 (50 m separation) 

Pulse Index 

Time index of detected peaks for each multipath 

Detected 
path #1 

Detected 
path #2 

Detected 
path #3 

Detected 
path #4 

Detected 
path #5 

1 

1115 

1349 

1405 

0 

0 

2 

1115 

1357 

1406 

1450 

1523 

3 

1115 

1353 

1402 

0 

0 

4 

1116 

1358 

1401 

0 

1526 

5 

1116 

1355 

1408 

1451 

0 

6 

1116 

1351 

1401 

0 

1529 

7 

1116 

1350 

1408 

1451 

1524 


Time index of detected peaks for each multipath 
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Pulse Index 

Detected 
path #1 

Detected 
path #2 

Detected 
path #3 

Detected 
path #4 

Detected 
path #5 

8 

1116 

0 

1408 

1451 

1524 

9 

1116 

1352 

1408 

1451 

0 

10 

1116 

0 

1413 

1443 

1524 

11 

1116 

1358 

1408 

0 

0 

12 

1116 

1351 

1410 

0 

0 

13 

1116 

1351 

1403 

1445 

1524 

14 

1116 

1352 

0 

1452 

1519 

15 

1117 

1354 

1404 

0 

0 

16 

1117 

0 

0 

1453 

0 

17 

1117 

1354 

1412 

1453 

0 

18 

1117 

1355 

1406 

0 

1524 

19 

1117 

1352 

0 

1453 

0 

20 

1117 

1352 

1405 

0 

1526 

21 

1118 

1355 

0 

1451 

1527 

22 

1118 

1352 

1409 

1449 

0 

23 

1118 

1352 

0 

1453 

1527 

24 

1118 

1352 

1412 

1453 

0 

25 

1118 

1352 

1405 

1454 

1525 

26 

1118 

0 

1412 

1448 

1529 

27 

1118 

1352 

1409 

1450 

1519 

28 

1118 

1352 

1414 

0 

1529 

29 

1118 

1357 

0 

0 

0 

Correlation 
Coefficient ( E j ) 

1 

0.9999987 

0.9999977 

0.9999978 

0.9999980 

Skew (F ( ) 

1.0000450 

0.9999948 

1.0000869 

1.0000336 

1.0000205 

Skew (Fpj ) 

-45 

5 

-87 

-34 

-21 

Skew error 

±3 

±23 

±32 

±35 

±36 
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Detected Peaks of 2 ms LFM chirp rom experiment #3 (50m) 
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Figure 32. The peak times for an experiment in which clock skew is present: The peak 
times u jt derived from the ensemble of n = 30 measurement signals depicted in 

Figure 30 are plotted in this figure. Note that each of the m = 6 multipath arrivals 

a.j , j = 1, 2,.. m has a corresponding ensemble U j = jw^. j ,j = l,2,...,m 

associated with it. 


The gradual decrease in stability of each successive detected path is also 
illustrated in Tables 4 and 5. Unlike the first path which is most likely the direct path 
between the two nodes, the following ones bounced off either or both the surface and the 
bottom of the lake, inducing variability in the arrivals. 

Figure 33 represents the received matched-filtered 10 ms CW signal recorded in 
experiment #3. This demonstrates the importance of selecting a probe signal which can 
be easily distinguished by the receiver. Although there is a general trend in the shape of 
each response, it is difficult to identify a constant peak over all of them. 
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First 15, 10 CW beacons (experiment #3) 



Figure 33. First fifteen 10 ms long CW beacon at 50 m (experiment #2). 

Finally, Figure 34 shows the impact of selecting a longer LFM pulse. In this case 
the main first peak is still clearly identifiable between each chirp, but its general shape is 
larger and the following chirps are less distinguishable. The combination of the minimal 
delay difference between the arrivals and the relatively large pulse length probably 
creates too much superposition to clearly recover subsequent arrivals after the first one. 
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First 15, 50ms long LFM chirps (experiment #3) 



Time sample index 2000 15 


Figure 34. First 15, 50 ms long LFM chirps at 50 m (experiment #3) 

Overall, the results from the experiment demonstrate time variability in the tested 
channel. The first arrival is always present and strong. However, the following arrivals 
are observed to be gradually less stable as they increasingly interacted with the channel 
boundaries. The algorithm is able to identify various multipaths based on the detected 
arrivals and determine the most stable one based on its correlation coefficient. The time 
indices of these arrivals clearly illustrate the presence of a skew between the two clocks 
as the time index of the main arrival slowly shifts to the right, as shown in Tables 4 and 5 
and the associated Figures 31 and 32. The variance of the results in Table 3 could be due 
the short period of time over which the chirps are transmitted. 
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VII. CONCLUSION AND RECOMMENDATIONS 


In this research we addressed the problem of clock synchronization in an 
underwater acoustic network. This is a very actively researched area for every type of 
communication channel. This is particularly challenging in the underwater acoustic 
environment due to its high latency characteristics, and it requires the development of a 
specific clock-synchronization protocol. Most of the approaches found in the literature 
assumed constant propagation delays between two nodes and ignore time variability of 
the channel. An overview of how and why time variance exists in the underwater 
acoustic environment was given and a method to evaluate and mitigate its effects as part 
of a clock-synchronization protocol was proposed. 

The proposed method is divided in two stages. The first stage evaluates the 
frequency drift between the clocks of two adjacent nodes (called the “skew”), and the 
second stage estimates the offset between the same two nodes. 

The first stage was simulated using a computer model and tested in a small, 
shallow, fresh water lake. Both simulated and experimental results show the 
effectiveness of the proposed method. However, further testing needs to be done on the 
basis of the following recommendations to address the variability in the estimated skew 
calculated from the experimental results. 

The second stage could only be simulated since an experiment would require 
modifying the firmware of the modem, which was not readily possible with the Teledyne 
Benthos acoustic modems. The results from the simulations are in accordance with the 
predictions from the theoretical model of the proposed method for offset correction. 

Additional experiments are required to fully develop the method and select the 
optimal parameters. Once completely tested, the addition of this capability to an 
underwater acoustic network will keep the clocks of the nodes synchronized. This 
synchronization should improve the utility of the network when deployed or used for 
long periods of time. 
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The following recommendations do not necessarily result from the simulation and 
experiments. Some of them are simply aspects that were omitted due to time or resource 
constraints. 

A. REPRODUCE EXPERIMENT FROM LITERATURE 

Most of the literature pertaining to underwater clock synchronization does not 
describe physical experiments in an actual underwater environment. The simulations 
reported in [1] and [7] should be reproduced using the protocol explained in this thesis. 
This would allow for some comparison among the various protocols and confirm that 
accounting for time variability improves the synchronization of the clocks. 

Another approach to compare the different protocols is to reproduce our 
experiment using the protocols from [1] and [7]. 

B. PSEUDORANDOM BINARY SEQUENCE (PRBS) 

Another type of pulse used in channel sounding is the pseudorandom binary 
sequence (PRBS) [10]. Just as the LFM chirp, the PRBS is characterized by its narrow 
ambiguity function, making it an excellent choice for measuring the impulse response of 
a channel. Its structure allows it to be more tolerant to superposition. The generation of a 
PRBS signal is not complicated and reproducing the experiment from Chapter VI would 
require minimal effort. 

C. LONGER RPI 

As discussed in Chapter VI, the variability of the results shown in Table 3 might 
originate from a short RPI. The experiment should be reproduced using a constant 2 ms 
LFM chirp while increasing the RPI duration. This would confirm our hypothesis that the 
results would improve with a longer RPI. 

D. OFFSET CORRECTION EXPERIMENT 

The simulation proves the theory of the offset correction algorithm. However, an 
experiment in the underwater environment can possibly highlight any neglected factors as 
was shown when testing our skew-correction algorithms. As stated earlier, technical and 
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time limitations prevented us from testing our offset correction algorithm in Del Monte 
Lake. An experiment should be designed and conducted to test the developed algorithm 
for offset correction described in Chapter IV and simulated in Chapter V. 

E. MORE REALISTIC ENVIRONMENT 

Del Monte Lake was selected mostly for its accessibility. As stated in Chapter III, 
the effects of time variability are location dependent. Since a shallow fresh water lake is 
not a probable environment where an underwater acoustic network would be deployed, 
the experiment should be reproduced in an environment more typical of where a time- 
synchronization protocol would be required and beneficial. 

One of the advantages of Del Monte Lake was its low noise level. This will not 
necessarily be true in an ocean environment. Future work should involve folding 
measurement noise into the model for skew and offset. Also, noise reduction algorithms 
should be explored. 

F. DETECTED ARRIVALS SORTING ALGORITHM 

The current algorithm is vulnerable to noise and superimposed arrivals. Future 
work should explore new algorithms to cope with these issues. 
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